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Abstract. Streamer discharges pose basic problems in plasma physics, as they are 
very transient, far from equilibrium and have high ionization density gradients; they 
appear in diverse areas of science and technology. The present paper focuses on the 
derivation of a high order fluid model for streamers. Using momentum transfer theory, 
the fluid equations are obtained as velocity moments of the Boltzmann equation; they 
are closed in the local mean energy approximation and coupled to the Poisson equation 
for the space charge generated electric field. The high order tensor in the energy flux 
equation is approximated by the product of two lower order moments to close the 
system. The average collision frequencies for momentum and energy transfer in elastic 
and inelastic collisions for electrons in molecular nitrogen are calculated from a multi 
term Boltzmann equation solution. We then discuss, in particular, (1) the correct 
implementation of transport data in streamer models; (2) the accuracy of the two term 
approximation for solving Boltzmann's equation in the context of streamer studies; and 
(3) the evaluation of the mean- energy- dependent collision rates for electrons required 
as an input in the high order fluid model. In the second paper in this sequence, we 
will discuss the solutions of the high order fluid model for streamers, based on model 
and input data derived in the present paper. 

PACS numbers: 52.25.Dy, 52.25.-y, 52.65.Kj, 52.25.Fi, 52.25. Jm 

Submitted to: J. Phys. D: Appl. Phys. 
1. Introduction 

Streamers are rapidly growing filaments of weakly ionized plasma, whose dynamics are 
controlled by highly localized space charge regions and steep plasma density gradients. 
The dynamics of the streamer ionization front are governed by electron dynamics 
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in electric fields above the break-down value; therefore the plasma is very far from 
equilibrium, and the neutral gas stays cold while the ionization front passes. 

Streamers occur widely in pulsed discharges, both in nature and in technology. 
As their size scales with inverse gas density, streamers occur in the limited volumes of 
plasma technological devices mostly at high gas densities, while planetary atmospheres 
also can host huge discharges at very low gas densities; so-called sprite discharges [H [21 13] 
are streamers, they exist at altitudes of up to 90 km in our atmosphere, i.e., at 
pressures down to well below 10 /ibar. For a cluster issue on streamers, sprites 
and lightning, we refer to [4J and the 19 original articles therein, discussing common 
issues of streamer dynamics in atmospheric discharges and plasma technology, including 
discharge evolution and subsequent chemical reactions. Streamers are a key element in 
gas processing in so-called corona reactors [HIE] as well as in high voltage technology [7j; 
they are used for the treatment of polluted gases [8] and water [9] or for plasma enhanced 
vapor deposition [10]. The plasma bullets observed in plasma jets in noble gases have 
now been identified with streamers [H], [12], [131 [lH [13 [IE]. Studies of breakdown 
phenomena at atmospheric pressure in short non-uniform air gaps [T7] or in microwave 
fields [IH] benefit from related streamer studies. Streamers appear as well in resistive 
plate chambers used in high-energy physics, where the streamer mode of operation must 
be carefully controlled and optimized [HI [201 [21] . 

The progress of streamer simulations with fluid models has recently been reviewed 
in [22]. The main point for the present paper is that essentially all these numerical 
studies are based on the classical fluid model, as we have called it in [231 [211 125] . 
The classical fluid model for the electron density contains a reaction term for impact 
ionization and attachment, a drift term accounting for electron displacement in the local 
field, and electron diffusion; the reaction term is typically taken as a function of the local 
electric field. The structure of this classical model with reaction, drift and diffusion is 
based on basic symmetry considerations and conservation laws; it has emerged as a 
purely phenomenological model in the course of the past century. 

On the other hand, cross-sections for the collisions of electrons with molecules have 
become available with growing accuracy [211 [23 [23 [23 [30] . They are the input either for 
the Boltzmann equation or for Monte Carlo models. These microscopic particle models 
have not yet been appropriately linked with the much older phenomenological fluid 
models. This becomes evident when comparing their solutions. A negative streamer 
ionization front in nitrogen represented by a Monte Carlo model cannot appropriately 
be described by a classical fluid model even if the transport and reaction coefficients for 
the fluid model are derived from the Monte Carlo model [31]. In [31], the coefficients for 
the fluid model were derived as so-called bulk coefficients; therefore in this model the 
electron swarms evolve correctly, and ionization fronts in a given electric field propagate 
with the correct velocity, as they are so-called puled fronts [31]. However, the electron 
density behind the ionization front is too low. As we will see in section [3] of the present 
paper, a classical fluid model evaluated with flux coefficients gives better ionization 
densities behind the front, but too low front velocities. 
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A phenomenological extension of the classical fluid model with a gradient expansion 
considerably improves the approximation of the microscopic electron dynamics [321 133] , 
and it cures the deficiency of the local field approximation by including a dependence 
of the impact ionization rate on the local electron density gradient [23]. That the 
extended fluid model matches the Monte Carlo model well up to the moment of 
streamer branching, can also be seen in a recent comparison of three-dimensional 
model simulations in [25]. On the other hand, the model is still based on a local 
field approximation, therefore electron energies in the streamer interior erroneously are 
modeled as relaxing instantaneously to the low field in the streamer interior. 

A full simulation by Monte Carlo models would deliver the physically most reliable 
results, but this is computationally extremely costly. As a compromise between accuracy 
and efficiency, hybrid methods [25l [23l [Ml [34] have been developed in the past years to 
track the fast and energetic processes in the ionization front with a particle model, and 
to model the many electrons at lower fields in the streamer interior with a fluid model. 
Fluid models are also the most efficient for the long time evolution in the streamer 
interior where slow chemical reactions and a slow thermalization of the deposited energy 
sets in. 

These demands on fast and quantitative streamer simulations ask for the 
development of a fluid model that accurately incorporates the microscopic electron 
dynamics contained in Boltzmann or Monte Carlo models. We will here follow the 
strategy to derive such a model not from phenomenological considerations, but through 
a systematic derivation from the Boltzmann equation. 

This paper is the first in an ongoing investigation of high order fluid models for 
streamer discharges; it is focused on the derivation of first and higher order fluid 
models and on the derivation and correct implementation of transport data in these 
models. Section [2] describes the derivation of models of different order. The starting 
point is the set of balance equations obtained as velocity moments of the Boltzmann 
equation. Then the derivation of the classical first order model is briefly described, 
and previous approaches to higher order models are summarized. In section 12.41 our 
new high order fluid model is derived in a systematic manner. The balance/moment 
equations are closed after the balance equation for the energy flux. This is done by 
approximating the high order tensors in the energy flux balance equation by a product 
of two lower moments while the collision transfer terms are evaluated using momentum 
transfer theory. In section 3, electron transport properties as an input in fluid models 
are calculated using a multi term theory for solving the Boltzmann equation [351 ISS] . 
We pay particular attention to the accuracy of the calculation and the proper use in 
both the first and high order fluid models. The results of our multi term solution 
of Boltzmann's equation are compared with those obtained by the publicly available 
Boltzmann solver BOLSIG-I- developed by Hagelaar and Pitchford [37] for electrons 
in molecular nitrogen. BOLSIG-I- is a popular Boltzmann solver based on a classical 
two term theory so we have been motivated to check its accuracy and integrity against 
the advanced and highly sophisticated multi term Boltzmann solver developed by the 



High order fluid model for streamer discharges 



4 



group from the James Cook University and their associates ESI ESI EH]- In section 
4 we give the results for negative planar fronts obtained with the first order model with 
particular emphasis upon the consistent implementation of transport data. The results 
for various streamer properties obtained with different type of input data are compared. 
A thorough analysis of numerical streamer front solutions with our high order model, 
and a comparison with Monte Carlo results is contained in an accompanying second 
paper [76] . 

2. Derivation of first and higher order fluid models 

2.1. General considerations: Boltzmann equation and moment equations 

Our starting point is the Boltzmann equation for charged particles in an electric field E 



Here fi{r, c, t) is the distribution function in phase space at the position r and velocity 
c for each charged component i, V is the differential operator with respect to space r 
and Vc with respect to velocity c, Cj and mi are charge and mass of species z, and t is 
time. The right-hand side of equation ([1]), J{fi,fo), describes the collisions of charged 
particles with neutral molecules, accounting for elastic, inelastic, and non- conservative 
(e.g. ionizing or attaching) collisions, and /o is the velocity distribution function of the 
neutral gas (usually taken to be Maxwellian at temperature Tg). 

Streamer discharges have a characteristic nonlinear coupling between densities of 
charged particles and electric field. The space charge modifies the field, and the field 
determines the drift, diffusion and rate coefficients. The electric field has to be calculated 
self-consistently with Poisson's equation 



where eo is the dielectric constant, and and Ui are the charges and densities of species i 
that can be electrons and ions. As ions are much heavier than electrons, their motion is 
typically neglected within the ionization front, and the analysis focuses on the electron 
evolution. Thus, in what follows we suppress the charged particle index i in equation 
(II]) and focus in the electron dynamics. 

After solving equations ([1]) and ([2]), quantities of physical interest could be obtained 
as velocity 'moments' of the distribution functions, starting with the number density 
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followed by higher order quantities 
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furnishing the average velocity v 



(c), average energy 
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Equations ([T]) and represent a non-linear coupled system of partial differential 
equations in 6- dimensional phase space and in time, with a complicated collision 
operator J. For a complete problem description, appropriate initial and boundary 
conditions on /(r, c,t) in phase space must be implemented. 

For streamers, it is very difficult to solve equations ([H) and ([2]) due to the large 
gradients of electric fields and of charged particle densities at their fronts. Moreover, 
the streamers create a self-consistent field enhancement at their tips which allow them 
to penetrate into regions where the background field is too low for efficient ionization 
processes to take place. From a kinetic theory point of view, this is a non-stationary, 
non-hydrodynamic and non-linear problem where space and time should be treated 
on an equal footing with velocity in equation ([1]). Clearly, the numerical solution of 
equations ([1]) and ([2]) for the full streamer problem is a formidable task and a fluid 
equation treatment is more tractable. In a fluid approach, the problem of solving the 
Boltzmann equation for / in phase-space is replaced by a set of low order approximate 
(velocity) moment equations of / [351 SOI SI] ■ 

Fluid equations may be derived either directly as moments of ([1]) or from 
first principles using physical arguments. Following the first approach, the set of 
moment/balance equations can be found by multiplying ([1]) by 0(c) and integrating 
over all velocities 

dt (n(0(c))) + V ■ (n(c0(c))) - n^E ■ (Vc0(c)) = , (5) 

where <> represents the average over the velocity c of the charged particles, and is 
the collision term in the balance equation: 

C4> = -J <P{c)J{f)dc. (6) 

To derive the term (Vc0(c)) with its negative sign, a partial integration over c has been 
performed. 

If one now takes 0(c) consecutively equal to 1, mc, |mc^ and |mc^c, etc., equation 
(|5]) generates an infinite series of equations, a full solution of which would be equivalent 
to calculating / itself. In practice, one truncates after a certain moment equation; in 
this process obviously some information of the Boltzmann equation is lost. Therefore we 
now discuss the derivation of first and second order fluid models with their truncations 
of the moment equations, and then we derive systematically our new high order fluid 
model. 



2.2. The first order or classical fluid model 

First we derive the classical fluid model as an approximation from the Boltzmann 
equation. We use continuity and momentum balance equation, i.e., 0(c) = 1 and mc, 
and truncate the set ([5]) at the momentum balance equation, 

dtn + V- (nv) = Ci , (7) 
dt{nmv) + V ■ {nm{cc)) — neE = Cmc , (8) 
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where v = (c) is the average local electron velocity. Now the velocities are decomposed 
into an average velocity v plus random velocities c — v with zero mean {{c — v) = 0). 
On introducing the pressure tensor 

P = nm{{c — v)(c — v)) , (9) 
equation becomes 

dt (nmv) + V ■ {nmvv) + V ■ P — neE = Cmc (10) 
where the following identity was used 

V ■ {nm{cc)) = V ■ [nm{vv + v{c — v) -\- (c — v)v + ((c — v) {c — v)))] 

= V ■ {nm vv) + V ■ {nm ((c — v) {c — v)))] (11) 
The second term in the left-hand side of ( |T0|) can be expanded as 

V ■ (nm vv) = nm (t) ■ V) f + f V ■ (nmv) . (12) 

Now we substitute ( |T2l) into ( ITOl) . use the continuity equation and introduce the 
convective time derivative 

! = « + V, (13) 

which measures the rate of change in a reference frame moving with the mean drift 
velocity v of the electrons. The momentum balance equation then obtains the form 
dv 

nm— = neE -V ■ P + Cmc - rrivCi. (14) 

The physical interpretation of this equation is straightforward: the rate of change of 
the mean electron velocity is due to the force of the electric field plus forces due to 
the pressure of the electron fluid itself and due to internal forces associated with the 
coUisional interactions with a large number of neutral gas molecules at rest. It should 
be emphasized, however, that the form with the convective derivative d/dt is not useful 
for the analysis of a full streamer problem where local fields and therefore local mean 
electron velocities v vary largely in space and time. 

Because the system has been truncated, the yet unspecified tensor P (|9]) of electron 
pressure appears on the right-hand side of (HM . and a closure assumption needs to be 
found. If the distribution of random velocities is close to isotropic, the diagonal terms 
of P are equal and given by the scalar kinetic pressure p, 

P^pI = nkT, (15) 

where I is the unity tensor, k is the Boltzmann constant and T is the temperature. 
It should be emphasized here that the isotropy of the velocity fluctuations and of the 
pressure is a strong assumption at the streamer tip where the electric fields are high and 
strong pressure gradients exist. In this streamer region the random spread of electrons 
along the field direction can differ significantly from the perpendicular direction. 

The next simplifying assumption concerns the collision terms. An expression often 
used for momentum transfer by collision is 

Cmc = -nmuesv , (16) 
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which assumes that the force per unit volume exerted on the electrons due to 
collisions with neutral molecules is proportional to the average electron velocity. The 
proportionality constant is called the effective momentum transfer collision frequency; 
it accounts for momentum exchange in elastic and inelastic collisions. With this 
simplifying assumption and neglecting the transfer of momentum in non-conservative 
coUision^jl relative to other momentum transfer collisions (which is usually a good 
approximation), the momentum balance equation (|T^ becomes 
dv 

nm— = ncE — Vp — nmvesv. (17) 
dt 

If the rate of momentum change {dtv) /v is smaller than the rate of momentum transfer 
z/gfr, and if the gradients in electron energy can be neglected, one gets the following 
expression for the average flux of the electrons 

T = nv = ni2E - DVn, (18) 

where mobility and diffusion constant are given by 

/x=^, D = ^, (19) 

mZ/eff mZ/eff 

if the system is close to equilibrium. In this case the Nernst-Townsend- Einstein relation 
D kT 

- = — 20 

jj, e 

is valid. The steady-state form of equation f lT7|) is an acceptable approximation because 
the effective time constant for momentum transfer ^ at atmospheric pressure is 
generally much less than the time scale on which the local electric field varies within a 
streamer [311 112] . 

Further from equilibrium as in the head of the streamer where the electric fields are 
high, the approximation f|T5l) is not valid as the velocity fluctuations of the electrons are 
clearly anisotropic. Let us consider once more the special case where the average velocity 
of electrons is stationary. From equation (IT^ and neglecting momentum transfer in non- 
conservative collisions, we obtain 

mnuesv = enE — VP, (21) 

and as a result for the flux we have 

T = nfiE - —V ■ [n{{c - v) {c - v))]. (22) 

Comparing the last equation and equation ( ITSl) . it is obvious that quantity 

„^ ((c — v) (c — v)) , , 

D"" = ^ ^ (23) 

is the reminisce of the diffusion tensor often used in the drift-diffusion approximation 
instead of the diffusion constant. However, the physical interpretation of this quantity 
is not a priori clear. Though this quantity assumes the anisotropic nature of 

I Non-conservative collisions do not conserve particle numbers as they account for ionization, 
attachment or recombination reactions. 



High order fluid model for streamer discharges 



8 



the temperature tensor, it reduces to the diffusion constant D when the effective 
momentum transfer coUision frequency z/gg is independent of the electron energy. For an 
energy-dependent effective momentum transfer colhsion frequency, the straightforward 
generahzation of the Einstein relation for the diffusion constant (fT9|) to the diffusion 
tensor D or application of D* in the equation for the average flux (|T8l) is misleading 
and reader is referred to [l3l HU HS] where the so-called generalized Einstein relation 
were introduced using the momentum and energy balance equations. However, the 
generalization of the diffusion constant D appearing in ( |T8i) to diffusion tensor D is 
a very welcome step in fluid modeling of streamer discharges due to the often strong 
anisotropic nature of the diffusion tensor for certain gases and due to deviations of the 
electron velocity distribution function from a Maxwellian distribution of velocities in 
different spatial regions of the streamers. 

When the above approximation for the particle flux is inserted into equation ([7]), 
we get the well known reaction drift diffusion expression for the charged particle motion 
in the discharge 

dtn + V ■{^^{E)En- D{E)-Vn) = Ci. (24) 

Originally this approximation was not derived from the Boltzmann equation, but 
derived on purely phenomenological grounds. It is clear that the equation must have 
the structure of a continuity equation with a source term for electron generation and 
loss. The parametric dependence of the source term is open at this point; within the 
classical model the source term is assumed to depend on local particle densities and on 
the local electric field. The second approximation concerns the electron velocity v. On a 
time scale much larger than the collision frequency, the electron motion is overdamped, 
so the electrons must drift against the direction of the electric field according to the 
second term in ( 12^ . The stochastic and undirected motion is modeled in an ad hoc 
manner through the diffusion term; the fact that the diffusion correction is added to the 
drift term and not included in any other functional manner is not a priori clear and can 
actually be deduced from the above analysis. 

In conclusion, the lowest level of fluid approximation, also called classical model or 
first order model, is given by the equations 
dfh 

— = V-{D-Vn) + V- ifinE) + n{uj - ua) , (25) 

coupled to the Poisson equation for the electric field, 

= — («-«P -^n), E = -V(l). (28) 
Co 

Here Up and n„ are positive and negative ion densities, while z// and z/^ are the ionization 
and attachment collision frequencies due to electron-molecule collisions, and (p is the 
electric potential. The numerical solution of the system requires the transport 
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properties /x, D, vj and va as a function of the local electric field for the gas in question. 
The derivation and implementation of transport data is discussed in section 3. The 
numerical solutions of planar streamer fronts with these transport data are discussed in 
section 4. 

2.3. Second order models including the energy balance equation 

We now turn to models that include the second velocity moment of the Boltzmann 
equation, i.e., the energy balance equation. We will see that the closure of this equation 
is not a straight forward process, and we will discuss some approximations made in the 
literature. The results of a second order model and of our high order model for planar 
streamer fronts will be compared in our second paper [76] . 

The first important steps beyond first order fluid models of streamer discharges, 
to our knowledge, were carried out by Abbas and Bayle [16| HTj, and by Bayle and 
Carenbois |1H]. They employed a second order model which involves the energy balance 
equation to explore the zone at the streamer tip where the electron energy is not 
determined anymore by the local electric field only. However, as pointed out by Kanzari 
et al. [12], the accuracy of their model was restricted by the drastic assumptions taken 
in their energy balance equation. They evaluated the source term in the energy balance 
equation and corresponding averages by assuming a Maxwellian distribution for the 
electrons. Guo and Wu [19] have developed a more sophisticated second order model 
in which the Langevin theory was used to simplify the collision source terms with a 
priori knowledge of the relaxation times of electron energy and momentum. Kanzari 
et al. [12] have made an important step further by calculating the source term in 
the energy equation in a more consistent way where the ad hoc assumptions for the 
distribution function were avoided. A similar approach was later used by Eichwald et 
al. [50] in their simulations of the streamer dynamics and of the radical formation in a 
pulsed corona discharge used for flue gas control. Their results showed that the average 
electron velocity in the streamer head is largely overestimated by the classical first order 
model. As a consequence, electron density and radical density in the ionized channel 
were up to 50% higher than with the second order model. The salient feature of their 
theory is the fact that the heat flux term in the energy balance equation was explicitly 
neglected, but this, unfortunately, is of questionable accuracy. We will illustrate in 
streamer simulations in our second paper [76] that the energy flux plays an important 
role in determining the streamer behavior. Thus, particular care should be taken with 
the closure through specification of the energy flux. Even for charged particle swarms 
under non-hydro dynamic conditions one must be careful how to specify the heat flux 
vector [35], HOI HH [51]. We now discuss how the fluid equations should be closed for 
streamers, while we stress that the method itself is applicable to a much wider range of 
phenomena. 
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2.4- High order fluid model 

We now derive a fluid model including the energy flux equation, i.e., we truncate at the 
moment equation with the third power of velocity. We will argue that the energy flux 
equation is crucial for the success of the fluid model for streamers. Inserting different 
moments of c into equation ([5]), one finds 

_ + V-nt; = Ci, (29) 

d 

— {nniv) + V ■ {nra{cc)) — neE = Cmc , (30) 

d f /I 2\A , w ^ /I 2 



at \^\2 / J + ^ ■ (^^(2"''' V J - ^ei; ■ z; = , (31) 

^ (^n(^mc2c)^+V-(^n(^mc2cc)^-ne£;-^^(imc2c)^ = Ci^,.,.(32) 

Different forms of these equations that some readers might be familiar with, are given 
in the appendix. It will be shown in section 12.4.31 that the explicit form of the high 
order tensors appearing in the divergence term of the energy flux equation ( 132|) is not 
required. 

This system of equations is exact at this stage, but not very useful and the reason 
is twofold. First, there are more unknowns than equations, the familiar problem of 
closure. The set of fluid equations can be increased to an arbitrary size merely by 
taking higher velocity moments of the Boltzmann equation. Second, the definition of 
the system requires the respective collision terms Cs. 



2. 4-1- Collision processes. In contrast to Monte Carlo simulations or kinetic equations 
in which cross sections for charged particle scattering enter into the calculations in 
a fairly clear way, in fluid equations the collisions are treated in a variety of ways 
which are not necessarily consistent with the data itself and/or with the system under 
consideration. There are many examples which illustrate this issue and the reader 
is referred to [351 SOI E] for a detailed discussion. In the context of streamer and 
breakdown studies beyond first order fluid models, it has become common practice 
to evaluate the collision terms and averages by assuming some particular form of the 
velocity distribution function, usually Maxwellian jUJ HOl |52l [531 El] • 

Let us consider the form of the collision terms in the balance equations (I29ll-(l32l). 
In this work we characterize the elastic, inelastic and non- conservative collisions by 
corresponding average collision frequencies. A collision frequency z/(fr) for collisions 
between charged particles and gas molecules is related to the cross section a{vr) 
characterizing the process by 

Z/(t>r) = ngVr<j{Vr) , (33) 



High order fluid model for streamer discharges 



11 



where is the relative velocity and rig is the number density of the neutral gas. Further, 
we deal with weakly ionized systems where the interactions of charged particles with 
one another and with excited molecules is negligible. 

Our calculation takes into account that the charged particles can loose energy and 
momentum even in elastic collisions with the gas molecules as the finite mass and the 
thermal energy of the molecules are taken into account. The momentum exchange in 
inelastic collisions is included and z/^ = Vmivr) denotes the total momentum transfer 
collision frequency while Vi and z/^^ are inelastic and superelastic collision frequencies 
for the inelastic channel i. The total collision frequencies for attachment and ionization 
are denoted by va = ^A^Vr) and ui = viivr), respectively. We consider only single 
ionization with ionization energy e/, but the resulting ion can be left in any one of 
its internal excited states, characterized by an excitation energy Ae^ ^ and a collision 
frequency uf^. 

2.4-2. Momentum transfer theory. The momentum transfer theory has a long history 
in kinetic theory of gases and has proved very successful for describing charged 
particle transport in gases under the influence of electric and magnetic fields. It is 
discussed comprehensively in the textbook of Mason and McDaniel [Hj and in references 
[35| HOl SH HHl EHl EEl ET] . The main ideas can be briefly summarized as follows: 
(1) We need to express the average collision frequencies as a function of the mean energy 
in the center-of-mass reference frame. Thus, we replace the variables 



where e is the energy measured in the center-of-mass reference frame and /x^ is the 
reduced mass. 

(2) If we assume that the dominant contributions to the averages come from regions 
near the mean energy e, and that //(e) varies sufficiently slowly with e, then a Taylor 
expansion 



is expected to be a reasonable approximation. For conservative collisional processes such 
as elastic and inelastic scattering, only the first term of the expansion (1361) is considered. 
However, when energy- dependent non-conservative processes such as ionization and 
electron attachment are operative then the derivative term in (1361) becomes the leading 
term and must be included. 

(3) Using the momentum transfer approximation, the balance equations (I29l) - (l32l) get 
the form 



in expressions for collisional frequencies 



(34) 




(35) 




(36) 




(37) 



High order fluid model for streamer discharges 12 
d 

— {mnv) + V ■ {mn{cc)) — neE = —firnUmV — mnv[i'i + ^i^^] , (38) 
ot 

- - vn^. - neV^ - n ^ Vf^ , (39) 

i i 

^(ri^) + V- (^(^rnc^ccj^ - neE - (^^{\^^^^)^ = -nUm^ , (40) 

where e is the average electron energy, ^ is the average electron energy flux and mo is 
the mass of gas molecules. ( is given by |l5l 156] : 

2 fl_ / 2\ 1 I |2 



3 m + mo \ / / / 

The average collision frequencies for momentum and energy transfer 

Umie) = ngJ—am{e) , (42) 

V /^r 

De{e) = -^^u.^{e) , (43) 
m + mo 

are prescribed functions of the mean energy in the centre-off-mass frame 

^ = — ^ , 44 

m + mo 

where k is Boltzmann's constant and To is the neutral gas temperature. Other 
collision frequencies describing the inelastic, non- conservative and superelastic collisions 
appearing in equations (137|) - (H0|) are also functions of e. It should be emphasized that 
the equations of continuity (!37|) . of momentum balance (138|1 and of energy balance (139!) 
are valid for charged particles of arbitrary mass while the energy flux equation ( l40ll is 
obtained in the approximation of m/mo ^ 1. High-order corrections in momentum 
transfer theory corrections (e.g., high-order terms in the Taylor expansion ( l36l) ) could 
be added on the right-hand side of the system (137I) - (140|) if desired, without in any 
way changing the generality of the physical arguments associated with the closure 
assumptions presented below. 

2.4-3. Solution regimes and closure assumptions. The closure of the system of 
equations (I37l)- (l40l) requires approximations or assumptions on the form of the pressure 
tensor. The standard approximation for light particles such as electrons is that the 
pressure tensor can be taken as a scalar at the fluid level of approximation [35 | HO t BT | H5] . 
This means that, e ~ e ^ |mf ^, and the pressure tensor simplifies to 

2 

P = nkT^-neI, (45) 
3 

where T is the so-called temperature tensor that characterizes energy fluctuations 
even if the system is not in thermal equilibrium. This form of the pressure tensor 
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was employed in all previous swarm oriented studies [3S1 SOI SIl SSI EB E51 ES], as 
well as in the recent fluid models of streamer discharges ^2} [50] . However, if the 
velocity distribution signiflcantly deviates from isotropy in velocity space, then this 
approximation is problematic. For ions, the distribution function in velocity space is 
always anisotropic (even if elastic collisions between ions and molecules are predominant) 
and hence any assumption on an isotropic pressure tensor is wrong for ions. This problem 
can be avoided by considering the temperature tensor balance equation but this in turn 
contains further unknowns. The reader is referred to [57] for how to treat charged 
particle swarms under spatially-homogeneous hydrodynamic conditions. In streamer 
studies, the ions are usually considered as immobile or they are modeled by a reaction- 
drift-diffusion and local fleld approximation which is a reasonable approximation on the 
time scale on which a streamer ionization front passes a given point in space. 

The next step in the closure of the system of equations ( l37l) -( l40l) concerns the energy 
flux balance equation (HOl) . The third term can be simplifled as follows: 

2 \ 9c 1 ndC /.^N 

— — -mc c = mc— c H — mc — — = mcc H — mc 1 , (4o ) 

oc\2 / oc 2 oc 2 



where / is the unity tensor. Assuming again as in (H5l) that the temperature tensor is 
isotropic, and hence that 

(cc) ^ Y /, (47) 
and after some algebra, the energy flux equation can be written as 

^«) + V ■ (^n(^^mc^cc^^ - ^E(^^ne^ = -nu^^ . (48) 

The second term in the energy flux balance equation (HS!) or (HO!) is the divergence of 
the fourth power of the velocity averaged over the velocity distribution, (c^cc), while all 
other terms in the equation contain only the third power of velocity. Therefore this term, 
that we will call the quartic tensor, requires either the next equation in the sequence 
of moment equations or a closure approximation. This closure assumption must be 
physically transparent and consistent with the general structure of the equations. One 
way is to approximate the relevant term by a product of lower moments as 

c'cc)^p{c'){cc)^P^^e'l, (49) 

where we used f HTl) for the second equality. /3 is a parametrization factor, generically 
close to unity, when the higher order correlation term (c^cc) — (c^) (cc) can be neglected. 

With these closure assumptions the system of fluid equations ( j37j) - (l40l) for the 
electrons becomes: 
Ori 

— + V-nv = -n{uA - J^i) , (50) 
|(n^) + Av(n.) - n^E = -nv[u^ + + , (51) 
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d 3 
— {ne) + V ■ «) - cE ■ {nv) = -nu^ U - -kTo 



n 



^^""^^ + ^('^S^') " ^""^'^ = "^"^^''^^ • ^^^^ 
The parameter (3 that appears in the energy flux equation is a quantity close to unity, but 
can be used to fit the neglected higher order correlations. Variations of this parameter 
will be discussed in our second paper [76] . 

The system of equations fl5(l - flS^ has the following properties: (i) If the parameter 
/3 is specified and if the collision terms are given then the system of equations contains 
no further unknowns can be numerically solved for the density, average velocity, average 
energy and energy flux, (ii) The equations contain mean-energy collisional rates that 
should be carefully derived and implemented as elaborated in the next section, (iii) 
Attachment enters the momentum and energy balance equation in terms of the derivative 
of the attachment collision frequency while for ionization only the ionization collision 
frequency is present, (iv) The above equations are set for a single component gas; 
the generalization to gas mixtures proceeds through the generalization of the collision 
term to a sum of terms appropriately weighted according to the mole fractions of the 
respective species. 



3. Transport and reaction data 

3.1. Evaluation of the transport data: Boltzmann equation analysis 

Here we discuss how to evaluate and implement electron transport properties in both 
first and high order fluid models. In plasma modeling, transport coefficients of electrons 
and/or ions may be obtained either from swarm experiments, from solutions of the 
Boltzmann equation, or from Monte Carlo simulations. We note that swarm transport 
data are usually given in the literature in the form of tables as a function of reduced 
electric field, E/uq [351110]. In fluid plasma modeling, however, several important issues 
have to be considered before directly implementing data from the literature. 

First, transport coefficients are defined far from boundaries, sources and sinks of 
charged particles where the so-called hydrodynamic conditions prevail [35l |36l [Ml [39l 
[581159]. Under hydrodynamic conditions, the space-time dependence of the distribution 
function is expressible in terms of linear functionals of n{r,t). A sufficient functional 
relationship between the distribution function f{r,c,t) and n{r,t) in the case of weak 
gradients is the well-known expansion 

oo 

/(r, c,t) = J2 f^'\^) (-V)%(r, t) , (54) 

s=0 

where f^'^\c) are tensors of rank s and denotes an s-fold scalar product. Direct 
application of transport properties measured/calculated in different experimental 
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arrangements where often non-hydro dynamic conditions are present exphcitly or 
imphcitly, is problematic and should be avoided. Typical example of swarm data 
obtained under non-hydro dynamic conditions are those measured/calculated under the 
steady-state Townsend (SST) conditions [601 EI]- Before direct application of SST 
data one must perform careful swarm analysis and convert the SST data into the 
hydrodynamic transport coefficients. Details of this procedure are presented in j62j . 

Second, care must be taken when non-conservative collisions are significant. In 
the presence of non-conservative collisions there are two sets of transport coefficients, 
the bulk and the flux [381 ES]. Assuming the functional relationship (1541) the flux F 
and source term Ci (usually denoted by S{r,t) in previous swarm oriented studies 
[351 l36l [581 159] ) appearing in the equation of continuity ([7]) can be expanded 

r(r, t) = nW^"^ - D^"^ ■ Vn (55) 
S{r, t) = 5(°)n(r, t) - 5^1) ■ Vn{r, t) + 5^^) : VVn(r, t), (56) 

where W^*-*^ and D*-*^ define, respectively, the fiux drift velocity and flux diffusion tensor 
while S'^^^ are expansion coefficients of the source term. Substitution of expansions (l55l) 
and (l56l) into the continuity equation ([7]) yields the diffusion equation 
on 

— + W-n-D: VVn = -RaU (57) 

which defines the bulk transport coefficients 

Ra = -S^'\ (58) 
W = + 5(1)^ (59) 

D = D^*^ - 5(2), (60) 

where Ra is the loss rate, W is the bulk drift velocity and D is the bulk diffusion tensor. 

The basic difference between the bulk and flux transport coefficients should now be 
apparent. The bulk drift velocity is displacement of the mean position of the electron 
swarm and it characterizes the motion of the total ensemble of electrons. The presence of 
the electric fleld results in a spatial variation in the energy throughout the swarm. Under 
such conditions, the presence of non- conservative collisions (ionization/ attachment) 
may lead to a change in the position of the center-of-mass of the swarm. This effect on 
the bulk drift velocity is denoted by S^^\ On the other hand, the flux drift velocity W^*^ 
represents the rate of change of the position of the center-of-mass due to the electric 
fleld only and can be interpreted as the mean velocity of the electrons. Likewise the flux 
diffusion tensor D^*^ represents the rate of spreading of the swarm due to the electric 
fleld E and gradients in density Vn. The presence of non-conservative collisions may 
result in the variation of Vn throughout the swarm and a subsequent variation in the 
rate of change of the mean squared width of the swarm. Such effects are expressed 
by the second rank tensor S^'^\ The most appropriate procedure in plasma modeling 
would be to use the experimental swarm data (e.g. bulk values) for the analysis of 
the validity of the cross section and then to calculate the flux quantities which are 
necessary as input data in fluid modeling. More about duality of the hydrodynamic 
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transport coefficients and their implementation in fluid modeling can be found in the 
references [28l |35l EH SQl EH [59] . 

Third, while the first order fluid model requires electron transport data as a function 
of local reduced electric field, what really appears in high order fluid models are mean- 
energy-dependent coUisional rates. Since momentum transfer theory is used to determine 
the collision terms in the fluid equations, the most appropriate procedure would be 
to systematically reduce these fluid equations down to the swarm limit assuming the 
hydrodynamic regime. From this set of equations, one can then find relationships 
between coUisional transfer rates and the mean energy, in a self-consistent manner. This 
method was employed in previous works of Robson and co-workers (see for example [ID] 
and references therein), but we apply here a slightly different approach. Instead of 
using the so-called generalized Einstein relation to determine the mean energy from the 
transverse diffusion coefficient [iQl US] , the mean energy is directly calculated from the 
multi term solution of Boltzmann's equation. The correspondence between the mean 
energy and E/no is then used to find the correspondence between the mean energy and 
other relevant transport data. The momentum transfer collision rate is obtained from 

T^m = ^ , (61) 

where /i(e) is the electron mobility which is here a function of the mean energy. The 
energy-transfer collision frequency is calculated from equation (H3|) . Thus our procedure 
of determining the electron transport data as an input for high order fluid models is 
entirely consistent with the work of Boeuf and Pitchford [51] . 

The electron transport data employed in this work are calculated using a multi 
term theory for solving the Boltzmann equation. The methods and techniques are by 
now standard and the reader is referred to our previous works [351 EH [58]. Among many 
important aspects, we highlight the following important steps: 

• No assumptions on symmetries in velocity space are made, and the directional 
dependence of the phase-space distribution function in velocity space is represented 
in terms of a spherical harmonic expansion: 

oo I 

f{r,c,t) = Y,Y. fir,c,t)Y:^\c), (62) 

/=0 m=-l 

where Ym (c) are spherical harmonics and c represents the angles of c. In contrast to 
the classical two term theory (where the sum over / is performed only up to / = 1), 
the number of spherical harmonics is not restricted, and our method therefore is a 
truly multi term approach. The differences between two-term approximation and 
our full multi term approach for electron transport in nitrogen are illustrated in 
Section [31 The inadequacies of a Legendre polynomial expansion (when density 
gradients are not parallel to the field) are highlighted in our previous publications 
[351 [36] and avoided in this work. 
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• As discussed above, under hydrodynamic conditions a sufficient representation of 
the space dependence is an expansion in terms of powers of the density gradient 
operator. 

• The velocity (energy) dependence of the phase-space distribution function is 
represented by an expansion about a Maxwelhan at an arbitrary temperature in 
terms of Sonine polynomials. 

Using the appropriate orthogonality relations for the spherical harmonics and for 
the modified Sonine polynomials, the Boltzmann equation is converted into a hierarchy 
of coupled equations for the moments of the distribution function. These equations 
are solved numerically and all transport and rate coefficients are expressed in terms of 
moments of the distribution function [35l |36] . 

3.2. Cross sections and transport data 

In this section, transport and reaction coefficients for electrons in N2 at a temperature 
of 298 K are calculated as an input for ffist and high order fluid models. The ffist order 
model is based on the local field approximation; it requires mobility, diffusion coefficient 
and ionization rate as a function of the reduced electric field E/uq (where uq is the gas 
number density) . Compared to our previous work [58] , we extend the electric field range 
up to 3000 Td (1 Td=10~^^ Vm^). The high order fluid model requires average collision 
frequencies for momentum and energy transfer in elastic and inelastic collisions, and 
rate coefficients for all collision processes as a function of the mean electron energy. 

We use the cross sections for electron scattering in N2 provided by Stojanovic and 
Petrovic [66]. For elastic collision processes we use the original Boltzmann collision 
operator [67] while for inelastic processes we employ the generalization of Wang-Chang 
et al. [68]. The ionization collision operator is detailed in [39]. We assume that the 
ratio between the energy of the scattered electron and the total available energy in an 
ionizing collision is equally distributed between and 1; the same holds then, of course, 
for the ejected electron. We remark that at the high electron energies in the streamer 
tip, the assumptions on the energy division can considerably influence transport proffies. 
Furthermore, scattering is assumed to be isotropic. This can be problematic for high 
values of E/uq (generally for E/n^ > 1000 Td for electrons in N2) when electrons scatter 
predominantly in the forward direction [211 EH] • However, the errors in the calculated 
transport coefficients and rate coefficients are acceptable in fluid modeling for streamers 
in the range of the reduced electric fields E/uo considered in this work after appropriate 
renormalization of cross-sections for the scattering angle distribution [21]. 

We present both bulk and flux coefficients obtained by our multi term solution of the 
Boltzmann equation, and we compare them with the coefficients obtained by the public 
available Boltzmann solver BOLSIG+ derived from the same cross sections. BOLSIG-I- 
is based on the two term approximation p?] and provides exclusively flux transport 
data. 
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3. 3. Input data for the first order model 

Figure la shows the electron mobihty as a function of the reduced electric field E/uq. 
We observe that bulk and flux quantities start to differ visibly above a reduced field of 
approximately 150 Td, this means that the ionization processes start to be significant at 
this value of the field. As E/hq increases further, the effect becomes more pronounced, 
until the bulk mobility exceeds the flux mobility by approximately 30% at 3000 Td. The 
difference between bulk and flux mobility is the consequence of the spatial variation of 
the average electron energy within an electron swarm [SH [36l [39] . If the ionization rate 
is an increasing function of electron energy (as is the case for the parameters considered 
here), electrons are preferentially created in regions of higher energy resulting in a shift 
in the centre of mass position as well as in a modification of the spread about the centre 
of mass. In nitrogen up to 3000 Td, the electrons are preferentially created at the 
leading edge of an electron swarm and hence the bulk mobility is larger than the flux 
mobility. 

Figure la shows as well that the flux data obtained by our multi term solution of 
Boltzmann's equation and by the BOLSIG+ code agree well. Only for E/uq below about 
3 Td, the BOLSIG+ mobility is higher than our flux mobility. As we have successfully 
compared our multi term results with Monte Carlo results that include the thermal 
energy of the background molecules as well, the results shown in Fig. la suggest that 
the BOLSIG+ code should be carefully tested in the limit of thermal energies. 

Figure lb shows the diffusion coefficients as a function of the reduced field E/uq. 
The bulk and flux values of the longitudinal and transverse diffusion coefficients are 
compared with the isotropic diffusion coefficient calculated by the BOLSIG+ code. As 
for the mobility, flux and bulk data start to deviate for E/uq above approximately 150 
Td, indicating again the onset of ionization effects. The diffusion coefficients are more 
sensitive to the ionization processes than the mobility; the differences between bulk and 
flux data can reach almost 50% for E/uq approaching 3000 Td. For a more thorough 
analysis of the explicit influence of ionization processes on the diffusion coefficients one 
must consider the second order variations in the average energy along the swarm. This 
is beyond the scope of this work and we defer this to a future study. 

Figure lb clearly shows the anisotropy of the diffusion tensor, i.e. Dl ^ Dt', for 
the rage of E/uq considered here, the transverse diffusion coefficient is always larger 
than the longitudinal diffusion coefficient. This is due to the spatial variation of 
the average energy within the swarm and to the energy dependence of the collision 
frequency. The theory of anisotropic diffusion in an electric field is now textbook 
material (see for example [69]) and rather than a detailed discussion of the origin of 
this phenomenon, we prefer to highlight the implementation of the diffusion coefficients 
in fluid models of streamer discharges. For E/uq below approximately 30 Td, our 
results for the flux component of the transverse diffusion coefficient and those obtained 
by the BOLSIG+ code agree very well. For E/uq above 30 Td, however, BOLSIG+ 
with its two term approximation is clearly above our multi term solution. Here one 
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1000: , 





Figure 1. (a) Mobility, (b) longitudinal and transverse diffusion coefficient, and (c) 
ionization rate of electrons in N2 as a function of the reduced electric field E/uq, as an 
input for the first order model. Shown are flux and bulk data obtained by our multi 
term solution of the Boltzmann equation, and flux data obtained by the BOLSIG+ 
code. BOLSIG+ provides only an isotropic diffusion coefficient for panel (b). 



should bear in mind that BOLSIG+ treats the diffusion processes under the spatially 
homogeneous conditions and hence the transverse diffusion coefficient obtained under 
the spatially-inhomogeneous conditions should be used for comparisons. On the other 
hand, our one- dimensional fluid model requires the longitudinal diffusion coefficient as 
an input as we consider the spatial variations of the electron density and average electron 
energy only along the field direction. Therefore, particular care needs to be taken with 
implementation of the diffusion coefficients in fluid models of streamer discharges. 

Figure Ic shows that the ionization rate differs between the two term and the 
multi term calculations by up to 30 %. It is interesting to note that the two term 
approximation is less accurate in the energy region dominated by the vibrational 
excitation of N2 and for energies well above the ionization threshold. Surprisingly, for the 
electric field range between approximately 200 and 600 Td the two term approximation 
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increases in accuracy. In this energy region the cross sections for inelastic processes are 
much smaller than for elastic collisions. Similar but not identical observations have been 
made by Phelps and Pitchford |l70j . 

3.4- Input data for the high order model 




Figure 2. Average collision frequency (a) for momentum transfer v^m and (b) for 
energy transfer as a function of the mean energy of electrons in N2, as an input 
for the high order model. The three curves in each panel show the flux and bulk data 
obtained by our multi term solution for the Boltzmann equation and the flux data 
obtained by the BOLSIG+ code. 
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Figure 3. Rate coefhcients in N2 as a function of the electron energy, calculated 
with the multi term solution of the Boltzmann equations, (a) Rates for the cross 
sections listed in [BHl [70] and (b) momentum transfer rate, total inelastic rate and 
total ionization rate. 
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The average collision frequencies for momentum and energy transfer as a function 
of mean electron energy is shown in figure 2. Bulk and flux frequencies obtained with 
the multi term approach are compared flux values calculated with BOLSIG+. The flux 
values of the two approaches agree well. The differences between bulk and flux values 
of the average collision frequencies are a direct consequence of the differences between 
bulk and flux mobilities (see figure la). We see that the average collision frequency for 
momentum transfer is larger than for energy transfer. This follows from the fact that 
the energy loss in an elastic collision is proportional to the mass ratio m/mo ^ 1 (see 
equation (H3!) ) and thus much smaller than in an inelastic collision. 

Figure 3a shows the rate coefficients as a function of the mean electron energy 
for all collision processes. Figure 3b compares the momentum transfer rate, the total 
inelastic rate and the ionization rate. All results in figure 3 are derived with the multi 
term solution of Boltzmann's equation. 



4. First order streamer model with different transport data 

Here we present simulations of planar negative ionization fronts in N2 using the first 
order fluid model, while simulation results with the high order fluid model are presented 
in our second paper [76] . 



Numerical methods, initial and boundary conditions. 

In this subsection we briefly describe the numerical method, and the initial and boundary 
conditions used to solve equations fl25|) - fl28|) in one spatial dimension. The calculations 
are carried out in N2 at atmospheric pressure and at the ambient temperature of 298 
K. The ID simulations are started with the same initial Gaussian type distribution for 
electrons and ions 

(X - XqY 



n{x)\t=o = rii exp 



a2 



(63) 



in a gap parameterized with the coordinate x G [0,L], with L = 1.2 mm. We have 
chosen = 2 x 10^^ m"^, xq = 8 x 10^^ m and a = 2.9 x 10"^ m. The externally 
applied electric field is positive in the x direction, therefore electrons drift to the left. 
The field is fixed in the non-ionized region at the left boundary x = 0, providing a fixed 
electric field for the negative streamer ionization front to penetrate. In this work we 
consider reduced electric fields of 350 Td, 460 Td, 590 Td, 770 Td and 1000 Td. 

To investigate the sensitivity of streamer properties to the definition and accuracy 
of the transport data, we employ three different sets of data for /i, D and uj: bulk and 
flux data obtained by our Boltzmann equation analysis and the flux data obtained by 
the BOLSIG+ code. 

The finite volume method is used to spatially discretize the system (j25l)- (!28l) on a 
uniform grid with 1000 points. More details will be outlined in our second paper [76] 
where numerical methods for solving the high order fluid equations are presented in a 
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comprehensive way. To approximate the spatial derivatives in fl25p we use the second- 
order central difference discretization while the time derivatives are approximated with 
the Runge-Kutta 4 method. The continuity equation for the electron density has a 
second order spatial derivative, and therefore it requires two boundary conditions. For 
X = we use a homogeneous Neumann boundary condition {dxU = 0), so that electrons 
that arrive at this boundary may flow out of the system. For x = L we employ a 
homogeneous Dirichlet boundary condition {n = 0) to ensure that there is no outflow 
of electrons from the system. In any case, it should be noted that the electrons are well 
separated from the boundaries, therefore the actual boundary conditions do not matter. 
In equation ( |26ll the time derivative is approximated with the Runge-Kutta 4 method 

using the same time step as in equation fl25|l . In the ID case, equation fl28l) has the form 

g 

d^E = —{n - Uion) , (64) 

from which we can determine the electric field E by integrating over x and using the 
fixed value of at x = 0. 

4-2. Overview of simulation results with different transport data 

Figure H] displays the spatio-temporal evolution of electron and ion densities and of the 
electric field when the reduced electric field ahead of the front is fixed to 590 Td (or 
equivalently to 145 kV/cm in N2 at atmospheric pressure and temperature of 298 K). 
Calculations are performed for three different sets of transport data as indicated in the 
figures. We start with a Gaussian density distribution as described above. Though the 
transition from avalanche to streamer has been discussed many times within the past 
80 years [311 ED E2], the characteristics and main physical processes are discussed here 
to investigate the sensitivity to different transport data. 

In the early stage of evolution, we see that both the electron and the ion densities 
grow due to electron impact ionization. If this was the only mechanism, the electric field 
would remain unchanged and the ionization would continue indefinitely. However, the 
electrons drift in the direction opposite to the electric field while the positive ions would 
slowly drift in the opposite direction; as their mobility is so much smaller, this motion 
is actually neglected here and in most other streamer studies. As a consequence, the 
charge separation starts to distort the initially homogeneous electric field. Now figure H] 
shows that the ionization profiles at time 0.06 ns obtained with the bulk transport data 
are somewhat wider while their height is less than with our flux and BOLSIG+ data. 
This follows from the fact that the bulk mobility and diffusion constant are higher than 
the corresponding flux data and hence the centre of mass moves faster and the electron 
package spreads faster. As the ionization rate is the same in both cases, the height of 
the profiles obtained with the bulk data must be less than with the flux and BOLSIG+ 
data. As the evolution continues, the electric field in the ionised region gets completely 
screened, and further ionization processes cannot occur in this region anymore. The 
transition from avalanche to streamer is then completed. We mention in passing that 
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the complete screening of the interior field is due to the ID set-up and to the fact, that 
the field ahead of the front does not change in time. 
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Figure 4. Evolution of (a) tlie electron density, (b) the ion density and (c) the electric 
field in a negative planar ionization front in N2 with a reduced electric field of 590 Td 
ahead of the front. Shown are the spatial profiles obtained with three different sets of 
input data as indicated in the graph. 
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4-.3. Front velocities 

When the streamers have approached an approximately uniformly translating state, we 
see that the streamers obtained with bulk data propagate faster than those obtained 
with flux and BOLSIG+ data, in agreement with earlier studies [251 [311 [23] . section 
13.11 it was already discussed, that in general flux transport data should be used in 
fluid equations derived systematically from the Boltzmann equation; however, in the 
particular case of streamer ionization fronts with their pulled dynamics simulated with 
the half-phenomenological classical streamer model, the bulk coefficients approximate 
the front velocity better, with the drawback that the ionization level behind the front 
is too small [251 EH [23] . 

In figure [5] we display velocities of planar fronts as a function of E/hq. We compare 
the velocities obtained with different input data. The flux drift velocity as a function 
of E/uq is also shown. First, we see that the planar fronts move much faster than the 
electrons. This follows from the fact that the velocity of a planar front is the sum of 
the drift in the electric field at the front edge plus a term accounting for diffusion, for 
creation of additional electrons due to impact ionization and for the electron density 
profile [75]. The difference between front velocity and electron drift velocity increases 
with E/uq, up to more than a factor of 2 for the highest field displayed here. 

We remark that according to analytical theory [711 [75] . the front velocities in planar 
configurations (where the field does not decay ahead of the front) depend for a long time 
on initial conditions, and if the initial condition decays less than e"^*'"^'. A* = ^yuj/D, 
it will determine the front velocity for arbitrarily long times. For this reason, we here 
do not compare numerical with analytical results [75] . 

We finally note that the front velocities calculated with bulk transport data are up 
to 30 % higher than with flux data. This illustrates the sensitivity of the model to the 
input data. On the other hand, the velocities differ much less between our flux data and 
those obtained with BOLSIG+. To explore this issue in more detail, additional tests 
are required, particularly for atomic and molecular systems with large anisotropy of the 
velocity distribution function in velocity space. 

4.4- Ionization levels behind the front 

We now explore how the ionization degree behind the front depends on the transport 
data, and we compare with the analytical approximation 



This approximation becomes an identity, if diffusion can be neglected [73]; and it is 
an upper bound, if diffusion is taken into account [31]. This analytical result has been 
derived for planar fronts, and it is independent of the front velocity. 

In figure [6] we compare the ionization levels behind the fronts calculated with the 
three different sets of transport data and with the analytical upper bound, using our 
flux data. We see that the approximation (l65l) indeed serves as such a bound, but is 
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Figure 5. Velocities of planar fronts obtained by three different sets of input data 
as a function of the reduced electric field. The flux drift velocity of electrons is also 
included. 




Figure 6. Ionization levels behind planar fronts as a function of the reduced electric 
field, results are shown with three different sets of input data and with the analytical 
approximation (j65l) calculated with flux data. 



furthermore also a very good approximation of the numerical results, when using the 
same transport data. The errors of the two term approximation are negligible indicating 
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a weak sensitivity of the ionization level to the isotropy of the distribution function in 
velocity space that is assumed in BOLSIG+. The ionization level with bulk data is 
considerably lower. Generally, it is evident that the ionization level is much less sensitive 
to the type of transport data than the front velocity. Similar observations have been 
made in [3T] . 

5. Conclusion 

In this paper we have derived a high order fluid model for streamer discharges. Our goal 
has been to develop a comprehensive theory at a level of sophistication appropriate for 
the electron energy distributions very far from equilibrium and for the steep electron 
density gradients that characterize a streamer front. The first steps are reported 
in the present paper, which deals with theoretical foundations and phenomenology. 
After examining the application of the Boltzmann equation, we have proceeded to 
the evaluation of the velocity moments of Boltzmann's equation. To make contact 
with previous works and to place our theory in a much broader context, we have 
demonstrated how to derive the classical drift-diffusion approximation often used in the 
plasma modeling community to model streamer discharges. Then we have introduced 
our more sophisticated high order fluid approach, proceeding directly from Boltzmann's 
equation and systematically discussing the critical assumptions required to close the 
system of equations and to evaluate the collision terms involved. Momentum transfer 
theory has been used to approximate collision terms in the high order fluid model while 
high order tensors appearing in the energy flux equation have been specified in terms 
of previous moments. In contrast to previous works, it has been emphasized that the 
energy flux equation plays a pivotal role for the correct description of streamer dynamics. 
The fluid equations obtained as velocity moments of the Boltzmann equation have been 
closed in the local mean energy approximation and coupled to the Poisson equation to 
calculate the modification of the electric field by space charges. The numerical solutions 
of the high order fluid model for planar fronts and their discussion are deferred to the 
following paper. 

The second important aspect of this paper concerns the application of transport 
data in fluid models of streamer discharges. In order to illustrate this issue, we have used 
the first order fluid model to investigate the temporal evolution of negative planar fronts 
in pure nitrogen. We have focused on the way in which the inherent streamer properties 
such as the velocity of a streamer or the ionization level behind the front are influenced by 
different transport data employed as an input in fluid equations. Our primary goal was to 
show which aspects of kinetic theory developed for swarm physics and particularly which 
segments of data would be important for further improvement of streamer models. It was 
shown and illustrated that the direct application of transport data from the literature 
without knowledge of origin and nature of the data is problematic and can often lead 
to significant errors in the profiles of various streamer properties. In this respect, the 
origin and nature of transport data must be known and, if appropriate, suitably modified 
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before implementation in the fluid models. This is particularly important for collisional 
transfer rates required as input in high order fluid model. We have also discussed 
the validity of transport data obtained by a two term theory for solving the Boltzmann 
equation. Our general sentiment was that two term data may be used for fluid modeling 
of streamers only if demands for high accuracy are relaxed. 
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Appendix A. Alternative forms of the high order model 



To make contact with previous work [121 HTl HH, [50], we here present different forms 
of equations fl29|) - fl32|) . When we evaluate the averages over the velocities c in these 
equations and perform a considerable amount of algebra, we find 



dn „ ^ 
d 

— [nrav) + V ■ {nmvv) + V ■ P — neE = Cmc, 
at 
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(A.4) 



where P is the pressure tensor, q = ^nm{ (c — v)'^ (c — v)) is the heat flux vector, 
S = ^nm{{c — v)^{c — v){c — v)) is the high order pressure tensor and Q = 
nm{{c — v){c — v){c — v)) is the high order heat flux tensor. The pressure tensors 
P and S are the second order tensors while the high order energy flux tensor Q is a 
third order tensor. 

Using the continuity equation (lA.ip . the momentum balance equation (lA.2p can be 
written as 



nm 



V ■ P - enE = Cr, 



mvCi. 



(A.5) 
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This equation is equivalent to the equation ( IT^ and can be used to exclude the energy 
of macroscopic motion from the energy flA.3l) and energy flux (IA.4I) balance equations, 
respectively. For this purpose, we multiply momentum balance equations ( 1A.2P and 
(]A.5|) by V and after addition of one of the resulting equations to another one, we 
obtain 

^ {^nmv^^ + V ■ [Kimv'^v^ + (V ■ P) ■ v = t;C„e - (A.6) 

This is the balance equation for the energy of macroscopic motion. Taking this equation 

back into the energy and energy flux balance equations, we obtain an alternative and 
more compact form of the fluid equations 

^ + V-nv = C,, (A.7) 

at 

dv 

nm— + V ■ P - enE = C„vr, (A.8) 

Jtilp) + Ip^^ ■v) + V-q + P:Vv = (A.9) 
^ + q.Vv + q{V-v) + Q:Vv + V-S+[^-^E){T + ^pl) = (A.IO) 

where Vr = c — v is the random velocity. The pressure p is defined as one-third the 
trace of the pressure tensor 

P=IT.P^^ = J I i^-^ffdc, (A.ll) 

i 

while T is the stress tensor 

T = P-pI, (A.12) 

where / is the unity tensor (diagonal elements equal to unity). In the equations (1A.8I) - 
(lA.lOp . the convective time derivative f lT3|) has been utilized. The collisional terms on 
the right-hand-side of the balance equations (lA.8l) -( |Anol) are given by 

Cmvr = - / mVrJ{f)dVr, (A. 13) 



C^.m.^ = - j\rnvlJU)dVr^ (A.14) 

Cim^;2vr = - j ]^mvlVrJ{f)dVr, (A.15) 

where J{f) is the collision operator and / is the distribution function of charged 
particles. It should be noted that balance equations (IA.7P -f lArT0|) can be derived directly 
from the Boltzmann equation. In such case, however, the Boltzmann equation must be 
transformed into an equation for before the velocity moments are taken. 
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